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ABSTRACT 


In  this  thesis,  an  adaptive  two  dimensional  least  mean  squares  (LMS)  algorithm  and 
a  recursive  least  squares  (RLS)  algorithm  are  developed  from  the  one  dimensional  algo¬ 
rithms.  Design  of  the  two  dimensional  LMS  and  RLS  algorithms  are  studied  for  accu¬ 
racy  based  on  the  results  of  a  two  dimensional  system  identification  model  which  was 
used  in  testing  the  algorithms.  Application  of  the  two  algorithms  is  demonstrated 
through  computer  simulation  in  which  the  adaptive  filters  are  employed  in  a  noise 
canceler  and  an  adaptive  line  enhancer  and  applied  to  an  image  processing  problem. 
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I.  INTRODUCTION 

The  area  of  digital  signal  processing  has  experienced  a  rapid  growth  in  the  last  dec¬ 
ade.  Reasons  for  this  have  been  the  tremendous  advances  in  integrated  circuit  technol¬ 
ogy,  and  some  significant  developments  in  digital  processing  techniques  achieved  during 
this  period.  Included  in  these  developments  are  methods  which  extend  certain  one  di¬ 
mensional  digital  signal  processing  techniques  to  two  dimensions.  This  extension  is  by 
no  means  trivial.  Three  significant  factors  must  be  considered: 

1.  more  degrees  of  freedom  are  available  in  a  two  dimensional  system;  this  gives  a 
system  designer  more  flexibility  than  the  one  dimensional  case, 

2.  one  dimensional  problems  generally  involve  considerably  less  data  than  two  di¬ 
mensional  problems,  and 

3.  the  mathematical  methods  for  handling  two  dimensional  systems  are  generally  less 
complete  than  those  for  one  dimension. 

As  the  techniques  for  processing  multidimensional  data  have  improved,  the  appli¬ 
cations  of  digital  signal  processing  have  spread  from  one  dimensional  to  two  dimen¬ 
sional,  to  n-dimensional  data.  There  exists  many  physical  phenomena  that  inherently 
depend  on  two  or  more  independent  variables.  In  the  prediction  of  weather  and  in 
seismic  analysis,  the  data  depends  on  more  than  one  independent  variable.  Moreover, 
data  originating  from  one  dimensional  processes  can,  in  some  cases,  be  considered  two 
dimensional.  For  instance,  data  from  periodic  or  cyclic  processes  can  be  represented  as 
two  dimensional  arrays  by  using  their  periodicity. 

Besides  the  general  applicability  of  two  dimensional  signal  processing  in  the  above 
cases,  other  areas  which  have  experienced  significant  growth  in  recent  years  include  ra¬ 
dar,  sonar,  and  radio  astronomy.  The  two  dimensional  processing  of  images  is  also  a 
very  important  one.  Images  depend  on  two  spatial  variables,  and  are  continuous  in 
nature.  However,  if  we  digitize  them  and  assume  linear  models  in  their  formation,  dis¬ 
tortion,  and  recording  we  then  have  techniques  that  can  be  used  in  their  processing. 
Depending  on  the  applications,  different  processing  techniques  are  used,  notably:  en¬ 
hancement  and  restoration  of  images,  and  segmentation  and  encoding  of  images.  For 
details,  see  References  1,2,3,  4  . 

This  brief  discussion  about  applications  of  two  dimensional  digital  signal  processing 
shows  that  they  are  found  in  a  wide  variety  of  fields.  In  this  thesis,  we  are  interested  in 
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extending  one  dimensional  adaptive  filtering  techniques  to  two  dimensions  and  then 
applying  this  in  the  area  of  image  restoration.  Once  the  transition  from  one  to  two  di¬ 
mensions  is  understood,  the  extension  to  n-dimensional  signal  processing  is  fairly 
straightforward. 

A.  OBJECTIVES  OF  THE  THESIS 

The  use  of  the  Wiener  filter  has  proven  to  be  a  very  powerful  tool  in  the  area  of 
image  restoration  (Refs.  1,2,  4  J.  However,  one  disadvantage  is  the  fact  that  the  Wiener 
filter  operates  under  the  assumption  that  the  image  is  stationary  which  generally  is  not 
the  case  in  an  image  processing  problem.  In  one  dimensional  signal  processing  the 
Wiener  filter  concept  has  provided  a  basis  for  various  adaptive  filtering  algorithms. 
Within  these  adaptive  algorithms,  the  filter  possesses  charaaeristics  which  can  be  mod¬ 
ified  to  achieve  some  end  or  objective  and  is  usually  assumed  to  accomplish  this  "adap¬ 
tation”  automatically  without  the  need  for  substantial  intervention  by  the  user.  The 
adaptive  filter  can  "learn"  the  signal  characteristics  when  first  turned  on  and  thereafter 
can  track  changes  in  these  characteristics.  The  first  objective  of  this  thesis  is  to  examine 
a  two  dimensional  Wiener  model  for  image  restoration  which  can  then  be  extended  to 
adaptive  algorithms.  Due  to  its  relatively  low  computational  requirements  and  the  fact 
that  it  will  work  in  a  variety  of  signal  environments,  the  least  mean  square  (LMS)  algo¬ 
rithm  will  be  investigated  first. 

The  second  objective  of  this  thesis  is  to  develop  a  second  two  dimensional  adaptive 
algorithm.  In  this  case,  we  will  work  with  the  recursive  least  square  (RLS)  algorithm 
which  offers  faster  convergence  than  the  gradient-search-type  algorithms  but  generally 
involves  a  greater  cost  per  data  sample  and  more  numerical  difficulties. 

The  final  objective  is  to  implement  the  LMS  and  RLS  algorithms  within  a  system 
identification  model,  a  noise  canreler,  and  an  adaptive  line  enhancer.  Each  algorithm 
possesses  several  variants  and  various  parameters  w'hich  can  be  modified.  A  represen¬ 
tative  sample  of  the  various  outputs  will  be  exan^ed  and  the  results  compared  in  order 
to  see  which  provides  the  optimum  solution  under  given  conditions. 

B.  ORGANIZATION  OF  THE  THESIS 

Chapter  II  is  designed  to  review  the  development  of  the  one  dimensional  Wiener 
filter  and  then  extend  it  to  two  dimensions  where  it  can  be  incorporated  into  a  two  di¬ 
mensional  LMS  adaptive  algorithm.  Computer  simulation  of  the  algorithm  within  a 
system  identification  model  is  performed  and  the  results  are  shown.  The  two  dimen¬ 
sional  RLS  algorithm  is  derived  in  Chapter  111  and  again  the  results  of  the  algorithm 
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within  a  computer  simulated  system  identification  model  are  shown.  Chapter  IV  con¬ 
tains  the  results  of  implementing  the  LMS  and  RLS  algorithms  in  a  noise  canceler  and 
an  adaptive  line  enhancer.  Conclusions  concerning  the  results  are  also  presented. 
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II.  TWO  DIMENSIONAL  ADAPTIVE  LEAST  MEAN  SQUARE 

ALGORITHM 


In  this  chapter,  we  will  review  the  derivation  of  the  one  dimensional  Wiener  filter 
and  then  develop  an  algorithm  to  extend  it  to  the  two  dimensional  case.  We  will  use  the 
two  dimensional  Wiener  filter  to  achieve  a  two  dimensional  least  mean  square  (LMS) 
adaptive  filter  algorithm  which  will  be  demonstrated  to  be  useful  in  image  processing  by 
applying  the  algorithm  in  a  noise  canceler  mode  and  in  an  adaptive  line  enhancer  con¬ 
figuration  for  the  restoration  of  images  corrupted  by  noise. 

A.  ONE  DIMENSIONAL  WIENER  FILTER 

The  problem  of  estimating  one  signal  from  another  is  one  of  the  most  important  in 
signal  processing.  In  many  applications,  the  desired  signal  is  not  available  or  observable 
directly.  Instead,  the  observable  signal  is  a  degraded  or  distorted  version  of  the  original 
signal.  The  signal  estimation  problem  is  to  recover,  in  the  best  way  possible,  the  desired 
signal  from  its  degraded  replica.  One  typical  example  which  we  will  deal  with  in  this 
paper  is  an  image  recorded  by  an  imaging  system  that  has  been  corrupted  by  noise.  The 
problem  is  to  undo  the  noise  induced  distortion  and  restore  the  original  image. 

This  represents  the  classic  one  dimensional  problem  in  communication  theory  where 
we  must  obtain  an  estimate  of  a  signal  of  interest,  which  can  be  observed  in  the  presence 
of  some  additive  noise.  In  other  words,  the  available  information  about  the  signal, 
5(/i),  is  contained  in  the  received  signal; 

u{n)  =  s(m)  +  w{n)  (2.1) 

where  w{n)  is  the  noise.  Therefore,  we  must  process  this  available  signal  u{n)  through 
an  optimal  processor  that  produces  the  best  possible  estimate  of  s(rt). 

In  order  to  establish  a  two  dimensional  Wiener  filter,  we  must  first  develop  a  one 
dimensional  algorithm.  This  task  has  been  approached  from  many  different  directions 
[Refs.  5,6,7].  The  formulation  by  Haykin  [Ref  8]  provides  the  most  logical  extension  to 
two  dimensions.  First,  we  consider  a  tapped  delay  line  filter  similar  to  Figure  1  on  page 
5.  The  filter  consists  of  a  set  of  delay  elements,  a  corresponding  set  of  adjustable  tap 
gains  or  coefficients  /i(l),  h{.\f)  connected  to  the  tap  inputs,  and  a  set  of  adders 
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Figure  1.  Tapped  Delay  Line  Filter 

for  summing  the  resultant  outputs.  The  filter  is  driven  by  a  random  time  series  produc¬ 
ing  the  sequence  u{n),  u{n  —  1),...,  u{n  —  M  -I-  1)  as  the  M  tap  inputs  of  the  filter. 

We  can  express  the  signal  produced  at  the  filter  output,  >(n),  by  the  convolution 
sum: 


>•(«) 


M 


h{k)u{n-k^  1). 


(2.2) 


We  desire  a  filter  which  in  some  way  minimizes  the  difference  between  some  desired  re¬ 
sponse,  d{n),  and  the  corresponding  value  of  the  actual  filter  output.  This  difference  can 
be  denoted  as 


e{n)  =  d{n)  -  >(n)  (2.3) 

where  e{/i)  is  called  the  error  signal.  In  Wiener  theory,  the  filter  is  optimized  by  mini¬ 
mizing  the  mean-square  value  of  this  error  signal,  e{n)  . 
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Let  the  mean-square  value  of  the  error  be  denoted  by 


MSE=E{e\n)}  (2.4) 

where  £  {.}  is  the  expectation  operator.  This  mean-square  value  is  a  real  and  positive 
scalar  quantity,  representing  the  average  normalized  power  of  the  error  signal,  e{n). 
Substituting  Equation  (2.3)  into  (2.4)  yields 

MSE  =  £  {d\n)}  -  2  £  {^/(«M«)}  +  E  tv^(«)}-  (2-5) 

Next,  substituting  Equation  (2.2)  into  Equation  (2.5)  and  then  interchanging  the  orders 
of  summation  and  expectation  in  the  last  two  terms,  we  get 


M 


MSE  =  £  {</^(«)}  -  2  \h{k)  E  {din)u{n  -  /fe  +  1)) 


*1  msal 


h{k)h{m)  E  («(«  —  k-\-  !)«(«  —  m  +  1)}  . 


(2.6) 


Assuming  that  the  input  signal  u(n)  and  the  desired  response  d(n)  are  jointly  stationart', 
the  three  terms  on  the  right-hand  side  of  the  above  equation  may  be  interpreted  as  fol¬ 
lows: 


1.  The  expectation  £  {4’(m)}  is  equal  to  the  mean  square  value  of  the  desired  response 
d(n); 


P^=E{d\n)). 


(2.7) 


2.  The  expectation  £  {d(n)u(n  -  A  -f-  I)}  is  equal  to  the  cross-correlation  function  of 
the  desired  response  d(n)  and  the  input  signal  u(n)  for  the  lag  of  k-1.  We  can 
therefore  v.Tite  the  single  summation  term  on  the  righthand  side  of  Equation  (2.6) 
as  follows; 


hik)E{d{n)u{n  -  k  +  1)}  = 


h{k)p{k-l). 


(2.8) 


3.  Finally,  the  expectation  £  {u{n  —  k  +  l)u{n  —  m  -1- 1)}  is  equal  to  the 
autocorrelation  function  of  the  input  signal  u(n)  for  the  lag  of  m-k; 


r{m  —  k)  =  £  {«(n  —  k  +  1)u(m  —  m  +  1)}  .  (2.9) 


Accordingly,  we  can  rewrite  the  double  summation  term  on  the  righthand  side  of 
Equation  (2.6)  in  the  form 


^  h{k)h{m)  £{«(«  —  k  +  1)m(/i  —  m  -I- 1)} 
h[k)h{m)r{m  —  k) . 

csl  m»1 


ms] 


(2.10) 
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Thus,  substituting  Equations  (2.7),  (2.8),  and  (2.10)  into  Equation  (2.6),  we  find  that  the 
expression  for  the  mean  square  error  may  by  rewTitten  in  the  form 

M  M  M 

MSB  =  Pa  -2  ^hik)p{k  -  1)  +  ^  ]^A(A:)A(m)r(/n  -  k) .  (2. 1 1) 

By  differentiating  Equation  (2.11)  with  respect  to  h(k)  and  setting  it  equal  to  zero,  we 
have  the  following  set  of  M  simultaneous  equations 


p(/c-l) 


/io(m)r(m  -  k). 


1.2 . M. 


(2.12) 


This  system  of  M  simultaneous  equations  is  called  the  normal  equations  w'ith  optimum 
filter  coefficients  as  the  unknowns.  With  the  following  definitions. 


(2.13) 


(2.14) 


KO)  r(l) 

r(l)  r(0) 


r(M-l)  r(M-2) 


r{M-  1) 
r{M-2) 


riO) 


(2.15) 
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we  can  rewrite  the  normal  equations  of  Equation  (2.12)  in  matrix  form  as 


p  =  Rho.  (2.16) 

This  represents  the  discrete-time  version  of  the  Wiener-Hopf  equation. 

B.  TWO  DIMENSIONAL  WIENER  FILTER 

The  following  derivation  parallels  the  development  by  Hadhoud  and  Thomas,  how¬ 
ever  this  research  and  simulation  were  completed  separately  and  prior  to  the  publication 
of  Reference  9.  In  order  to  be  applicable  for  an  image  processing  problem,  we  must 
extend  the  formulation  in  the  previous  section  to  two  dimensions  [Refs.  10,11J.  This  is 
accomplished  by  developing  a  basic  two  dimensional  Wiener  filter  as  shown  in 
Figure  2  on  page  9.  Within  this  filter,  we  use  two  input  images,  the  reference  array  X 
and  the  primary  input  array  D.  The  primary  input  image  D  is  a  two  dimensional  array 
which  represents  the  ideal  image  plus  additive  noise,  while  the  reference  image  X  is  noise 
that  is  assumed  to  be  correlated  to  the  noise  in  the  primary  input.  Both  the  input  arrays 
are  MxM  in  dimension.  The  Wiener  filter  is  an  NxN  causal  FIR  filter  with  a  set  of 
weights  fVj  defined  as 


IF/0,0) 

IF/0,1)  ... 

IF/0,.V-1) 

IF/1,0) 

I^/l.l)  ... 

IF/l,rt-l) 

IF/ A' -1,0) 

IF/ /V- 1,1)  ... 

1 

1 

which  minimize  the  mean  of  the  squared  error,  e^,  between  the  filter  output  and  the  de¬ 
sired  input  D.  We  designate  J  as  the  iteration  number  given  by  J  =  mM  +  n  where  m  and 
n  take  on  the  values  from  0  to  M-1.  Just  as  in  the  one  dimensional  case,  the  filter  out¬ 
put,  y{m,n),  is  the  convolution  sum  of  the  filter  mask  and  the  reference  input  X  which  is 
given  by 


y{m,n)  = 


H'jm)X{m  -  l,n  -  k). 


During  the  jth  iteration  the  input  from  array  X  is  represented  by  A)  where 


(2.18) 
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Figure  2.  Two  Dimensional  Wiener  Filter 


AXm,n  —  1)  ...  — A'+l) 


(2.19) 


.Y(m  -  /Y  +  1,/j)  A'(m  -  +  I,n  -  I)  ...  ;ir(m  -  A^  +  I,n  -  N  +  1) 


Using  the  Equations  (2.17)  and  (2.19),  Equation  (2.18)  can  now  be  written  as 


i 


N-\  V-1 


(2.20) 


for  the  jth  iteration. 

This  output  y(fn,n)  can  now  be  subtracted  from  one  pixel  D(m,n)  of  the  array  D  to 
produce  the  error  signal  at  the  jth  iteration 


)V-l  ^1 


ej  =  D(m,«)  -  2  2  lVjil,k)Xim  -  l,n  -  k). 


(2.21) 


Since  the  purpose  of  the  Wiener  filter  is  to  provide  a  set  of  weights  which  minimize 
the  MSE,  we  can  denote  it  as 


MSE^Eief} 


(2.22) 


where  £{.}  is  the  expectation  operator. 

Using  a  mathematical  derivation  similar  to  the  one  dimensional  case  of  the  previous 
section,  we  can  see  that 


+ 


N-\ 

ej  =  D\m,n)  -  2^  Y  W'/Z.A)  Z)(m,rt)  J(m A) 

V-l  V-l  ,V-1  V-1  ^  ^ 

E  Z  Z  Z  -  P'^'  - 

fco  Ifco  «35 


(2.23) 


Substituting  Equation  (2.23)  into  Equation  (2.22)  yields 


^-1  iV-l 

MSE  =  EiD\m,ny\  -  2^  Ywj^Uk)  ElD{m,n)  X{m  -  l,n  -  ky] 

^=,1  ^=,1  ^1:^1 

+  Z  Z  Z  Y,^J^’^)^VJ<P^<l)E[.X{m-l,r,-k)X{,m-p,n-q)-]. 

fco  Iteo 


(2.24) 


Defining  P  as  the  crosscorrelation  matrix  between  the  desired  response  D{m,n)  and  the 
reference  input,  R  as  the  input  autocorrelation  matrix,  and  W  as  the  optimum  Wiener 
weight  matrix,  we  can  rewrite  Equation  (2.24)  as 


MSE  =  ElD\m,n)-\  -  2^  ^Wjil,k)  P{l,k) 

iV-l  A'-l  ,V-I 


y^Vjil,k)Wj{p,q)  R[p-l,q-k). 
Finally,  if  we  minimize  the  MSE  with  respect  to  WJ^l,k)  then  we  have 


N-\ 


qSi 


This  equation  is  the  two  dimensional  equivalent  to  Equation  (2.12). 
of  Equation  (2.26)  is  given  as 


P  =  RW 

where  the  elements  of  this  equation  are  defined  as  follows 


R  = 


W  W  [«2]  . 

m  W  . 

[^_,v+2] 

.  ™ 

w  = 

’  ^0  ' 

p  = 

'  Po  ' 
p. 

Pn-\ 

(2.25) 

(2.26) 

The  matrix  form 

(2.27) 

(2.28) 

(2.29) 

(2.30) 


II 


Within  the  R  matrix  in  equation  (2.28)  each  element  is  a  block  Toeplitz  matrix  repres¬ 
ented  by  the  equation 

R{p  -  l,q  -k)  =  EZX{m  ~  l,n  -  k)X{m  -  p,n  -  9)]  (2.3 1) 

where  l,k,p,q  range  from  0  to  N-1. 

C.  TWO  DIMENSIONAL  LEAST  MEAN  SQUARE 

One  means  of  obtaining  an  approximate  solution  for  the  optimum  weight  matrix, 
W,  is  the  use  of  a  two  dimensional  LMS  algorithm  which  is  depicted  in  Figure  3  on  page 
13.  This  differs  from  Figure  2  on  page  9  in  that  the  error  ej  is  used  to  update  the  filter 
coefficients  before  shifting  the  data  uindow  Xj  across  the  reference  input  for  the  next 
iteration.  As  in  the  one  dimensional  LMS  algorithm,  we  are  updating  the  coefficient 
matrix  by  adding  the  present  matrix  to  a  change  proportional  to  the  negative  gradient 
of  the  error  where  the  one  dimensional  instantaneous  estimates  of  the  gradient  vector 
are  based  on  sample  values  of  the  input  and  the  error  signal  e{n)u[n)  (Ref.  8).  For  the 
two  dimensional  jih  iteration,  we  define  the  updated  matrix  as 

=  (2.32) 

where 

updated  weight  matrix 
W/  previous  weight  matrix 

Pi  scaler  multiplier  controlling  the  rate  of  convergence  and  filter  stability 
{e){X^  estimate  for  the  2-D  instantaneous  gradient 

The  previous  equation  can  also  be  written  as 

=  \Vj^l,k)  +  2p  {ej)  X(m  -  l,n  -  k).  (2.33) 

These  two  equations  give  the  two  dimensional  weight  updating  algorithm  for  the  LMS 
adaptive  filter.  The  algorithm  we  have  developed  here  may  be  implemented  without  any 
form  of  matrix  operations,  averaging,  or  differentiation. 

The  value  of  p  may  be  chosen  based  on  the  desired  tracking  ability,  steady-state 
mean  square  error,  and  convergence  speed.  In  signal  processing,  there  are  several 
methods  for  determining  a  suitable  value,  however  in  many  of  these  cases  it  requires 
knowledge  of  the  eigenvalues  and  eigenvectors  of  the  autocorrleation  matrix.  In  image 


Figure  3.  T»o  Dimensional  Least  Mean  Square 


processing,  one  method  which  does  not  require  these  values  for  the  computation  of  n  is 
trial  and  error  based  on  output  image. 

An  alternate  method  used  in  one  dimensional  design  which  again  does  not  require 
a  priori  knowledge  of  the  autocorrelation  matrix,  discussed  by  Bitmead  and  Anderson 
[Ref.  12],  is  the  normalized  LMS.  Using  our  previous  notation,  this  method  can  be  ex¬ 
tended  to  two  dimensions  by  first  considering  the  two  dimensional  LMS  update  equation 
(2.32).  In  this  equation,  we  redefine  the  step-size  parameter,  fi,  for  a  given  /  and  k  as 


13 


nm  = 


(2.34) 


g 

P  +  a]{l,k) 


in  which  a,  represents  the  normalized  step  size  chosen  between  zero  and  two,  P  is  an¬ 
other  small  positive  constant,  and  aj{l,k)  is  one  of  the  values  from  the  matrix 


o]{0,0) 

Oj{0,l)  ... 

ffy(O.fV-l) 

o]H,0) 

0;"(l.iV-l) 

a]{N-\,0)  orj(N-U)  ... 

o]{N-\,N-\) 

(2.35) 


The  element  required  from  the  matrix  depends  upon  the  current  values  of  /  and  k  being 
used  for  equation  (2.33).  The  matrix  may  be  be  intitialized  with  the  values  of 
(X{m  —  /,«  —  ^))*  and  to  update  the  values  in  aj,  we  use  the  following  equation 

=  p{o]m)  +  {I-  p){X{m  -  l,n  -  k)f  (2.36) 


where  p  is  a  weighting  factor  between  zero  and  one.  This  ensures  that  the  value  of  p 
does  not  become  large  enough  to  cause  the  algorithm  to  become  unstable. 

D.  IMPLEMENTATION 
1.  System  Identification 

In  order  to  initially  test  the  tw’o  dimensional  least  mean  square  algorithm,  we 
established  a  system  identification  model  shown  in  Figure  4  on  page  15.  Within  this 
model  the  output  of  the  known  FIR  filter,  d{m,n),  is  defined  as 

d{m,n)  —  .4  W(m,n)  +  .6  W{m  —  1  ,n) 

-  .3  W{m,n  -  I)  +  .2  IF(m  -  1,«  -  1) . 


The  output  of  the  two  dimensional  least  mean  square  filter,  ,  consists  of  a  set  of 
adjustable  coefficients  and  is  defined  as 

y{m,n)  =  ^0  W{m,n)  +  Al  IF(m-  !,«) 

+  A2  W{m,n  —  1)  -t-  A'i  }V{m  —  l,n  —  1) . 


The  adaptive  filter  output  y(ni,n)  is  then  compared  with  the  known  system  output 
d(m,n)  to  produce  an  error  signal  e(m,n),  defined  as  the  difference  between  them. 
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KNOWN 

FILTER 


d(m , n) 


Figure  4.  T>vo  Dimensional  System  Identification 

The  operation  of  the  adaptive  filter  is  to  minimize  the  error  signal  by 

providing  an  adaptive  process  for  adjusting  the  coefficients.  For  this  adaptive  process 
we  use  tlie  update  equation  (2.32)  for  the  two  dimensional  least  mean  square  algorithm 
developed  in  the  previous  section 


(2.32) 

For  this  model  a  32x32  white  gaussian  noise  matrix  was  used  as  the  input  and 
the  rate  of  convergence  can  be  be  seen  in  Figure  5  on  page  16  and  Figure  6  on  page 
17.  Following  600  iterations  all  the  coefficients  had  converged  to  within  10*^  of  the  ac¬ 
tual  coefficient  value  and  the  error  was  3x10"’.  A  computer  program  for  this  system 
identification  model  is  given  in  Appendix  A. 

2.  Adaptive  Noise  Canceler 

The  usual  method  of  estimating  a  signal  corrupted  by  additive  noise  is  to  pass 
the  composite  signal  through  a  filter  that  tends  to  suppress  the  noise  while  leaving  the 


15 


signal  relatively  unchanged.  The  noise  canceler  and  adaptive  line  enhancer  developed 
by  Widrow  [Ref.  13)  arc  well-documented  ways  of  doing  that.  Adaptive  noise  canceling 
is  a  variation  of  optimal  filtering  that  is  highly  advantageous  in  many  applications.  It 
uses  an  auxiliary  or  reference  input  derived  from  one  or  more  sensors  located  at  points 
in  the  noise  field  where  the  signal  is  weak  or  undetectable.  This  input  is  filtered  and 
subtracted  from  a  primary  input  containing  both  signal  and  noise.  The  reference  input 
and  the  noise  in  the  primary  input  arc  therefore  correlated,  and  as  a  result  the  primary 
noise  is  attenuated  or  eliminated  by  cancellation. 
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Figure  6.  LMS  System  Ideiitiricatioii-Rate  of  Convergence 

The  basic  noise  canceling  concept  is  illustrated  in  Figure  7  on  page  18.  A  signal 
is  transmitted  over  a  channel  to  a  sensor  that  receives  the  signal  plus  noise,  w,.  The 
combined  signal  and  noise  s  +  forms  the  "primary  input"  to  the  canceler.  A  second 
sensor  receives  a  noise  n,  which  is  correlated  in  some  way  with  the  noise  n^.  This  sensor 
output  provides  the  "reference  input"  to  the  canceler  The  noise  n,  is  filtered  to  produce 
an  output  y  that  is  a  close  replica  of  ««  .  This  output  is  subtracted  from  the  primary 
input  s  +  Mo  to  produce  the  system  output  s  +  n^—y. 
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Figure  7.  Adaptive  Noise  Caiiceler 


In  the  sjsteni  shown  in  Figure  7,  the  reference  input  is  processed  by  our  two 
dimensional  least  mean  stjuare  filter.  Thus  the  filter  operates  under  changing  conditions 
and  will  readjust  itself  continuously  to  minimize  the  error  signal. 

1  he  computer  program  which  simulates  this  noise  canceling  model  is  provided 
in  Appendix  B.  We  will  discuss  the  results  and  conclusions  concerning  various  simu¬ 
lations  in  Chapter  4. 

3.  Adaptive  Line  Enlianeer 

A  special  case  of  adaptive  noise  canceling  is  when  there  is  only  one  signal  x(n) 
available  w'hich  is  contaminated  by  noise.  In  such  a  case,  the  signal  x{ft)  provides  its 
own  reference  signal  >•(«),  which  is  taken  to  be  a  delayed  replica  of 
x{n):y{/i)  =  x(n  —  A),  as  shown  in  Figure  8  on  page  19.  The  adaptive  filter  will  respond 
by  canceling  any  components  of  the  main  signal  x-(/j)  that  are  in  any  w'ay  correlated  with 
the  secondary’  signal  >■(/;)  =  x(n  —  A).  Suppose  the  signal  jf(n)  consists  of  two  compo¬ 
nents:  a  narrow'band  component  that  has  long-range  correlations  and  a  broadband 
component  which  w’ill  tend  to  have  short-range  correlations.  One  of  these  could  repre¬ 
sent  the  desired  signal  and  the  other  an  undcsired  interfering  noise.  Suppose  that  the 
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Figure  8.  Adaptive  Line  Enliaiicer 


delay  A  is  selected  so  that  it  falls  between  the  correlation  lengths.  Since  A  is  longer  than 
the  elVcctive  correlation  length  of  the  broadband  component,  the  delayed  replica  will  be 
entirely  uncorrelated  with  the  broadband  part  of  the  main  signal.  The  adaptive  filter 
will  not  be  able  to  respond  to  this  component.  On  the  other  hand,  since  A  is  shorter 
than  the  correlation  length  of  the  narrowband  component,  the  delayed  replica  that  ap¬ 
pears  in  the  secondary  input  will  be  correlated  with  narrowband  part  of  the  main  signal 
and  the  filter  will  respond  to  cancel  it.  Note  that  if  A  is  selected  to  be  longer  than  both 
correlation  lengths,  the  secondary  input  will  become  uncorrelated  with  the  primary  input 
and  the  adaptive  filter  will  turn  itself  off.  In  the  opposite  case,  when  the  delay  A  is  se¬ 
lected  to  be  less  than  both  correlation  lengths,  then  both  components  of  the  secondary 
signal  will  be  correlated  with  the  primary  signal,  and  therefore  the  adaptive  filter  will 
respond  to  cancel  the  primary’  x{n}  completely.  The  computational  algorithm  for  the 
adaptive  line  enhancer  is  shown  in  the  following  three  equations: 

vr  \f 

•^('0  =  -  "0  =  -  m  -  A)  (2.40) 

m=0  ni=0 
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ein)  =  xin)  -  x{n) 

h„{n  +  1)  =  h„{n)  +  2  n  e{n)x(n  -m- A)  m  =  0,1,2,...,A/ 


(2.41) 


(2.42) 

For  this  model  we  also  developed  a  computer  simulation  which  incorporates  our 
two  dimensional  LMS  algorithm  and  it  is  provided  in  Appendix  C.  The  results  and 
conclusions  will  again  be  discussed  in  Chapter  4. 
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III.  TWO  DIMENSIONAL  ADAPTIVE  RECURSIVE  LEAST  SQUARES 

ALGORITHM 


In  the  previous  chapter,  under  the  LMS  algorithm,  the  available  data  samples  were 
used  in  order  to  attempt  to  move  the  current  estimate  of  the  impulse  response  to  the 
optimum  value.  This  approach  has  the  advantage  of  being  simple  to  implement  but 
carries  with  it  the  disadvantages  that  it  can  be  slow  to  approach  the  optimal  weight 
vector  and,  once  close  to  it,  will  usually  fluctuate  around  the  optimal  vector  rather  than 
actually  converge  to  it  due  to  the  effects  of  approximations  made  in  the  estimate  of  the 
performance  function  gradient. 

To  overcome  these  difficulties,  we  examine  another  approach  in  this  chapter  which 
uses  the  input  data  in  such  a  way  as  to  ensure  optimality  at  each  step.  This  alternative 
algorithm  is  based  on  the  exact  minimization  of  the  least  square  criteria.  This  algorithm 
is  known  as  recursive  least  squares  (RLS). 

A.  ONE  DIMENSIONAL  RECURSIVE  LEAST  SQUARES 

As  in  the  LMS  case,  the  one  dimensional  RLS  algorithms  is  developed  in  several 
different  ways  (Refs.  5,  8  ].  Orfanidis  [Ref.  6)  provides  a  derivation  which  we  will  con¬ 
sider  prior  to  our  extension  to  two  dimensions.  The  tapped  delay  line  shown  in 
Figure  9  on  page  22  will  provide  the  reference  for  the  following  discussion.  We  begin 
by  replacing  the  LMS  estimation  criteria  of  MSE  —  £[«*]  by 


n 

-i 


MSE  =  ;  e%k)  k  =  0,1 . n 


where 


e{k)  =  x(k)  -  ^{k) 


and  x{k)  is  the  estimate  of  x{k)  produced  by  the  Mth-order  Wiener  filter 


X  =  2]^im)jik-m) 
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Figure  9.  One  Dimensional  RLS  Reference  (TDL) 


Substituting  equation  (3.2)  into  (3.1)  and  setting  the  derivative  with  respect  to  h  to  zero 
we  find  the  least-squares  analogs  of  the  orthogonality  equations 
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dMSE 

ah 


-2}e(k)y{k)  =  0 


which  may  be  rewritten  in  their  normal  equation  form  as  follows 


(3.5) 


(xik)  -  h^yikMk)  =  0 


y{k)y{k) 


Mk)yik) 


(3.6) 


(3.7) 


Defining  the  quantities 


n 


/?(«)  = 

Xy(^)y(^)'' 

(3.8) 

r(rt)  = 

n 

^^^)y(^) 

(3.9) 

we  can  then  write  the  normal  equation  as  R{n)  h  =  r(«),  with  solution  h  =  R{n)-'  r{n). 
Note  that  the  n-dependence  of  R{n)  and  r(n)  depend  on  the  current  time  n,  therefore, 

h(n)  =  R(/j)~'r(n)  =  R(n)r(n)  (3.10) 


where  P(n)  =  R(n)-'.  These  are  the  least  squares  analogs  of  the  ordinary  Wiener  sol¬ 
ution,  with  R(n)  and  r(«)  playing  the  role  of  the  covariance  matrix  £[y(n)y(«)T  and 
cross-correlation  vector  £[jr(«)y(rt)],  respectively.  The  RLS  algorithm  is  obtained  by 
UTiting  equation  (3.10)  recursively  in  n  and  then  using  the  following  matrix  inversion 
lemma  (Ref  5] 

(A  +  BCD)-^  =  A~'  -  A~'B{DA~'B  +  C“')"’Z>/<"’  (3.11) 

we  get  the  update  equation  for  the  P  matrix 

1  +  y{nfP{n  -  l)y(«) 
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Using  the  quantities  in  equations  (3.8)  and  (3.9)  to  satisfy  the  recursions  yields 


Ri»)  =  /?(«-!)  +  }’(«)>’(«)'■  (3.13) 

Tin)  =  r(«-I)  +  xin)yin)  (3.14) 

and  substituting  equations  (3.12)  and  (3.14)  into  (3.10)  after  some  mathematical  manip¬ 
ulations  we  find 


h(rt)  =  h(«  ~l)  +  Pin)  yin)  e(«)  (3.15) 

which  diflers  from  the  LMS  alogorithm  by  the  presence  of  the  P(n),  vice  fi,  in  front  of 
the  weight  correction  term.  Since  Pin)  =  /?(rt)“'  is  a  measure  of  the  covariance  matrix 
£ty(rt)y(«)  j,  the  presence  of  R(n)-^  makes  the  RLS  algorithm  behave  like  Newton's 
method,  and  hence  has  fast  convergence  properties. 

B.  TWO  DIMENSIONAL  RECURSIVE  LEAST  SQUARES 

Although  it  is  discussed  briefly  by  Wellstead  and  Caldas  Pintos  (Ref  14),  limited 
information  is  available  in  the  open  literature  in  this  area.  In  order  to  develop  a  two 
dimensional  recursive  least  square  (RLS)  algorithm  we  will  extend  the  one  dimensional 
adaptive  algorithm  developed  in  the  previous  section  in  a  method  similar  to  that  used  for 
the  LMS  algorithm. 

As  in  Chapter  2,  using  Figure  10  on  page  25  as  a  reference  we  see  that  the  basic 
filter  has  two  input  images.  The  primary  two  dimensional  input  array,  D,  is  the  ideal 
image  plus  additive  noise.  The  reference  image  X  is  the  noise  array.  Each  array  is  of 
dimension  M  by  M.  The  filter  mask,  W,  is  N  by  N. 

The  one  dimensional  RLS  algorithm,  equation  (3.15),  is  the  same  as  the  one  di¬ 
mensional  LMS  algorithm  with  the  exception  that  Pin)  replaces  /x.  Therefore,  if  a 
method  of  developing  a  P  matrix  for  the  two  dimensional  case  can  be  devised  we  can 
then  use  an  algorithm  similar  to  the  two  dimensional  LMS  algorithm  to  update  the  filter 
coefficients.  First,  let 
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Figure  10.  Two  Dimensional  RLS  Reference 
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Pj- 


P/0.0)  P/0.1) 


Pj{0,N^-l) 


P/,N^ -1,0)  Pj(N^-l,l)  P/^N^ -l,N^ -1) 


(3.16) 


which  represents  the  P  matrix  on  the  yth  iteration  where  j  is  defined  as  j  =  mM  +  n. 
As  in  the  two  dimensional  LMS  algorithm,  we  let  the  filter  coefficient  matrix  be  given 


as 


»'/0.0)  Wj{0,\) 

Wf\,0)  Wj,\,\) 


W^0,N-\) 


(3.17) 


H"/iV-l,0)  ir/A^-1,1) 


Wj^N-  l,iV-  1) 


and  the  input  data  window  as 


X{m,n)  X{m,n—\) 


X{m-N-\-  i,n)  X{m  -  N  +  \,n  -  1)  . 


Xim,n-N+  1) 


X{m-N+  \,n-N+  1) 


(3.18) 


In  order  to  establish  an  update  equation  for  the  P  matrix  in  which  each  element  will 
have  the  proper  dimensions,  we  use  equations  (3.17)  and  (3.18)  and  make  the  following 
transformations;  the  filter  mask,  equation  (3.17)  is  transformed  into  a  vector  defined 
as 
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and  we  transform  the  Xj  matrix,  equation  (3.18)  into 


XXj 


X(m,n) 
X(m,n  —  1) 

X{m,n  —  A''+  1) 


X{m  —  N  +  \,n) 
X{m  —  N  +  l,n  ~  1) 


(3.20) 


X{m-N+  \/t-N+  I) 


Using  the  quanitities  defined  above,  we  can  then  write  the  P  matrix  update  equation  as 


Pj.,{XXj){XXj')Pj_, 

\-^{XXj)Pj_,{XXj) 


(3.21) 


As  discussed  earlier,  the  one  dimensional  RLS  weight  updating  algorithm  difiers 
from  the  LMS  in  that  it  contains  the  P  matrix  while  the  LMS  algorithm  contains  n.  This 
must  be  considered  in  two  dimensions  since  P  is  a  matrix  and  /i  is  a  scalar.  Therefore 
using  the  previous  equations,  we  can  obtain  the  update  to  the  filter  coefficients  via 

=  WWj  -  {Pj){ej){XXj)  (3.22) 


where  is  given  as 


=  D(m,«)  -  y{m,n) 


(3.23) 


yim,n)  is  the  convolution  sum  of  the  image  in  the  reference  input  with  the  filter  window 


y{m,n)  = 


N-\ 


W^l,k)X^l,k) 


(3.24) 


and  the  value  Pj  is  updated  by  equation  (3.21). 

The  cost  function  of  the  RLS  criterion  could  be  modified  to  include  a  windowing 
of  the  input  data  as  follows 


MSB 


e\k) 


This  modification  results  in  the  following  modified  P  matrix  update: 


Pj 


W-i 


Pj_,{XXj){XXj)Pj_, 
\  +  {XXj)Pj_,{XXj) 


(3.25) 


(3.26) 


where  q,  the  averaging  or  "forgetting  factor"  is  a  positive  constant.  It  is  usually  chosen 
to  be  slightly  less  than  one,  thereby  diminishing  the  contribution  of  the  "older"  data. 
The  problem  of  data  nonstationarity  is  the  main  reason  for  introducing  this  type  of 
weighting  factor. [Ref.  6]  It  should  be  noted  that  the  usual  RLS  algorithm  is  attained 
when  q  =  1  and  that  no  changes  in  the  amount  of  computation  is  required. 


28 


C.  IMPLEMENTATION 


1.  System  Identification 

Following  the  mathematical  development  of  the  two  dimensional  recursive  least 
squares  algorithm,  we  implemented  and  tested  it  as  was  previously  done  with  the  two 
dimensional  least  mean  square  algorithm.  Our  initial  testing  was  done  in  a  system 
identification  model  as  shown  in  Figure  11.  As  in  the  two  dimensional  LMS  case,  we 
used  a  know'n  filter  with  the  following  output  equation 


d{m,n)  =  .41F(m,/i)  +  .61f(m—  !,«) 
—  1)  +  .2/F(/«—  l,n~  1) 


Figure  II.  Tito  Dimensional  RLS  System  Identification 


29 


The  two  dimensional  adaptive  recursive  least  square  filter  output  was  given  as 


y(m,n)  =  AO  lV{m,n)  +  A1  W{m  -  l,n) 

+  A2  IV(m,n  -  1)  +  A3  lV(rrr  -  1,«  -  1) 

The  error  signal  was  generated  by  taking  the  difierence  between  the  desired  (known) 
signal,  and  the  filter  output, y(m,w).  By  using  equations  (3.21)  and  (3.22)  to  up¬ 

date  the  P  matrix  and  the  filter  coefficients,  we  caused  the  filter  weights  to  move  toward 
their  optimum  values  and  the  error  signal  to  approach  zero. 

The  computer  program  for  this  model  is  provided  is  Appendix  D.  For  this  sys¬ 
tem,  a  32x32  white  gaussian  noise  matrix  was  generated  for  the  input  to  both  the  known 
filter  and  the  RLS  adaptive  filter.  The  rate  of  convergence  is  shown  in  Figure  12  on 
page  31  and  Figure  13  on  page  32  and  as  can  be  seen  after  approximately  70  iterations 
each  of  the  coefficients  had  converged  to  with  10-’  of  the  actual  value  and  the  error  was 
10-^ 

2.  Adaptive  Noise  Canceier/ Adaptive  Line  Enhancer 

These  two  systems  were  both  discussed  in  Chapter  2  with  regard  to  implemen¬ 
tation  of  the  two  dimensional  LMS  algorithm.  In  order  to  further  test  the  RLS  algo¬ 
rithm  we  also  implemented  it  within  the  noise  canceier  and  the  adaptive  line  enhancer. 
The  computer  programs  w'hich  were  used  for  the  simulation  are  given  Appendix  E  and 
Appendix  F,  respectively.  The  simulation  results  and  discussion  will  be  provided  in 
Chapter  4. 


Figure  12.  RLS  System  Identincatioii  Rate  of  Convergence 
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Figure  13.  RLS  System  Ideiitiricatioii  Rate  of  Convergence 
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IV.  RESULTS  AND  CONCLUSIONS 


In  this  chapter,  we  present  the  experimental  results  from  our  implementation  of  the 
two  dimensional  LMS  algorithm  derived  in  Chapter  II  and  the  two  dimensional  RLS 
algorithm  derived  in  Chapter  III.  The  two  algorithms  were  used  in  a  noise  canceler  and 
an  adaptive  line  enhancer,  as  discussed  in  Chapter  11,  to  solve  an  image  restoration 
problem. 

Figure  14  shows  the  original  image  of  a  house  which  is  comprised  of  128  by  128 
eight  bit  pixels.  The  image  has  a  mean  value  of  131.68  and  a  variance  of  3194.4.  In 
Figure  15  on  page  34,  we  have  corrupted  the  original  image  with  additive  white  gaussian 
.noise  (/’ero  mean  and  variance  1600). 

I - i 


Figure  14.  Original  Image 
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Figure  15.  Image  plus  Additive  Noise 


A.  NOISE  CANCELER  RESULTS 

In  order  to  solve  our  image  restoration  problem,  we  initially  implemented  both  the 
LMS  and  the  RLS  algorithms  in  a  noise  canceler  using  a  two  by  two  filter  matrix  and 
a  value  of  70x10  ’’  for  u  in  the  LMS  algorithm. 

Figu.''e  16  on  page  35  shows  the  output  when  the  LMS  algorithm  is  processed 
through  the  image  for  one  pass;  this  results  in  16,384  updates  of  the  filter  corresponding 
to  the  16,384  pixels.  Figure  17  on  page  36  is  the  results  following  two  passes  through 
the  image.  No  significant  improvement  was  realized  with  more  than  twm  passes. 

Implementing  a  two  by  two  RLS  adaptive  filter  in  the  noise  canceler  produced  very 
favorable  results  after  one  pass  through  the  image.  These  results  are  shown  in 
Figure  IS  on  page  36  and  compare  well  with  the  original  image  and  the  LMS  output 
after  two  passes. 

For  the  LMS  algorithm  our  best  results  were  acliieved  with  a  value  of  70x10  ’  for 
ii,  however  in  order  to  demonstrate  the  elTect  of  changing  this  value  Figure  19  on  page 
37  represents  the  output  with  /.<  -  35.rlO’^  and  Figure  20  on  page  37  is  produced  by  a 
spatially  varying  normalized  u  as  discussed  in  Chapter  IF  It  can  be  seen  that  fo,’'  \  a!'iou'< 


values  of  u,  dilTerent  features  within  the  image  were  restored  at  difFerent  levels.  In  Fig¬ 
ure  16  on  page  35,  the  edges  and  more  detailed  segments  of  the  image  were  restored 
better  than  in  Figure  19  on  page  37,  however  areas  of  similar  contrast  were  not  as  well 
restored.  When  using  the  normalized  ^  the  mean  value  more  closely  approached  the 
original  mean,  however  the  variance  was  reduced. 

Finally,  increasing  the  number  of  coefficients,  using  a  three  by  three  filter  matrix  in 
the  LMS  algorithm  we  obtained  the  results  shown  in  Figure  21  on  page  38  in  one  pass. 
This  filter  showed  no  significant  improvement  in  the  variance  compared  to  the  two  by- 
two  filter;  however,  it  did  improve  with  respect  to  the  mean. 


Table  I  on  page  38  shows  the  restoration  results  comparing  the  mean  and  variance 
of  the  various  outputs  with  that  of  the  original  image  and  the  image  plus  noise. 


Figure  16.  Restored  Image:  Noise  Canceler/LMS  (One  pas.s) 
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Figure  17.  Restored  Image:  Noise  Canceler/LMS  (Two  passes) 


Figure  18,  Restored  Image:  Noise  Canceier/RLS  (One  pass) 


Figure  19.  Restored  Image:  Noise  Canceler/LMS  (different  mu) 


Figure  20.  Restored  Image:  Noise  Canceler/LMS  (normalized  mu) 
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Figure  21.  Restored  Image:  Noise  Canceler/LMS  (3  by  3) 


Table  1.  NOISE  CANCELER  RESULTS 


TYPE  OF  FILTER 

MEAN 

VARIANCE 

Original  Image 

131.68 

3194.4 

Image  plus  Noise 

134.25 

1647.7 

One  pass  2x2  L.MS 

114.52 

2034.3 

Two  pass  2x2  LMS 

136.78 

3711.8 

One  pass  2x2  RLS 

139.04 

3607.7 

Onejjass  LMSIdifTerent 
mu) 

120.12 

1940.4 

One  pass 

LMS( normalized  mu) 

126.53 

1500.7 

One  pass  3x3  LMS 

119.13 

2017.3 
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B.  ADAPTIVE  LINE  ENHANCER  RESULTS 

We  next  implemented  our  two  algorithms  in  an  adaptive  line  enhancer.  In  order  to 
obtain  a  reference  input  from  the  primary  input,  we  used  a  two  dimensional  delay  oper¬ 
ator  of  (m-l.n).  As  in  the  noise  canceler  we  again  used  a  two  by  two  filter  matrix  and 
this  produced  the  results  shown  in  Figure  22  on  page  39  and  Figure  23  on  page  40  for 
the  one  pass  LMS  and  the  one  pass  RLS,  respectively.  No  significant  improvement  was 
noted  in  the  two  pass  LMS  compared  to  the  one  pass  LMS. 


Table  2  on  page  40  shows  the  comparison  of  each  output  mean  and  variance  \vjth 
the  mean  and  variance  of  the  original  image  and  the  image  plus  noise. 


Figure  22.  Restored  Image:  Adaptive  Line  Enhancer/LMS  (One  pass) 
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Figure  23.  Restored  Image;  Adaptive  Line  Enhancer/ RLS  (One  pass) 


Table  2.  ADAPTIVE  LINE  ENHANCER  RESULTS 


TYPE  OF  FII.TER 

MEAN 

VARIANCE 

Original  Image 

IM.68 

3194.4 

Image  plus  Noise 

134,25 

1647,7 

One  pass  2x2  LMS 

125.03 

2521.8 

Two  pass  2x2  LMS 

126,94 

2,529.9 

One  pass  2x2  RLS 

1 39.69 

2124.2 

C.  CONCLUSIONS 

In  this  thesis  we  have  extended  a  one  dimensional  LMS  and  a  one  dimensional  RLS 
adaptive  algorithm  to  two  dimensions.  As  we  examine  the  results  we  see  that  many  of 
the  one  dimensional  comparisons  between  LMS  and  RLS  are  applicable  to  two  dimen¬ 
sions: 

1.  the  rate  of  convergence  of  the  RLS  algorithm  is  in  general  faster  than  that  of  the 
LMS  algorithm  by  an  order  of  magnitude 

2.  this  superior  perfomtance  of  the  RLS  algorithm  is  attained  at  the  expense  of  a  large 
increase  in  the  computational  complexity 


3.  there  are  no  approximations  within  the  derivation  of  the  RLS  algorithm,  therefore 
as  the  number  of  iterations  approach  infinity,  the  least  squares  estimate  approaches 
the  optimum  Wiener  value 

4.  in  the  RLS  algorithm  the  correction  which  is  applied  is  computed  using  all  past 
data  whereas  the  LMS  algorithm  uses  only  the  instantaneous  sample  and  the  error 
signal;  this  is  not  necessarily  an  advantage  in  image  processing. 

The  objectives  of  this  thesis  were  successfully  completed.  Some  suggestions  for  fu¬ 
ture  work  include:  (i)  to  examine  these  two  dimensional  algorithms  in  other  applica¬ 
tions,  applying  them  to  areas  other  than  image  processing,  (ii)  to  extend  the  two 
dimensional  LMS  and  RLS  algorithms  to  n-dimensions  and  implement  them,  (iii)  to 
analyze  the  extension  of  these  one  dimensional  algorithms  to  two  dimensions  in  the  fre¬ 
quency  domain.  Multidimensional  digital  signal  processing  is  being  applied  in  many 
areas  today,  however  the  potential  for  future  growth  and  applications  appears  unlimited. 
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rj  o  I-*  no  n  o  o  n  o  n  o  o  o  o  o  o  o  n  ono 


APPENDIX  A.  LMS  SYSTEM  IDENTIFICATION  PROGRAM 

SYSTEM  IDENTIFICATION  OF  AN  ADAPTIVE  TOO  DIMENSIONAL 
LMS  FILTER  VS.  A  KNOWN  TOO  DIMENSIONAL  FILTER 

REAL  MU 

DIMENSION  A0(0: 32,0: 32),A1(0: 32,0:  32) ,A2(0:  32,0:  32) , 

*  A3(0:  32,0:  32) ,E(0:  32,0:  32) ,D(0:  32,0: 32),W(-1: 32,-1: 32), 

*  Y(0:32,0:32) 

******VARIABLE  DEFINITIONS****** 

A=N0ISE  AMPLITUDE 
S=N0ISE  VARIANCE 
AM=N0ISE  MEAN 
IX=SEED 

A0,A1,A2,A3=FILTER  COEFFICIENTS 
MU=SCALING  FACTOR  WITHIN  THE  LMS  ALGORITHM 
W=WHITE  GAUSSIAN  NOISE  GENERATED  BY  SUBROUTINE 
Y»FILTER  OUTPUT 
D=DESIRED  OUTPUT 
E=ERR0R(D  -  Y) 

*******initial  values******* 

A=l. 

MU=.  01 
A0(0,0)=0. 

A1(0,0)=0. 

A2(0,0)=0. 

A3(0,0)=0. 

S=l. 

AM=0. 

IX=65539 

*****0PEN  FILE  AND  FILL  INPUT  MATRIX  WITH  WHITE 
GAUSSIAN  noise******** 

OPEN  (UNIT=80,FILE='SYSID2  DATA’ ,STATUS=’ NEW’ ) 

DO  10  M=-l,31 
DO  20  N=-l,31 

IF((M.  LT.  0).  OR.  (N.  LT.  0))THEN 
W(M,N)=0. 

ELSE 

CALL  GAUSS(IX,S,AM,V) 

W(M,N)=V*A 
ENDIF 
CONTINUE 
CONTINUE 

******C0MPUTE  THE  UPDATED  FILTER  COEFFICIENTS 
C  USING  THE  TOO  DIMENSIONAL  LMS  ALGORITHM***** 

DO  40  K=0,31 
DO  50  J=0,31 

D(K,J)=. 4*W(K,J)+. 6*W(K-1,J)-.  3*W(K,J-1)+.  2*W(K-1.J-1) 
Y(K,J)=A0(K,J)*W(K,J)+A1(K,J)*W(K-1,J)+A2(K,J)*W(K,J-1)+ 


I 
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*  A3(K,J)*W(K-1,J-1) 

E(K,J)=D(K,J)-Y(K,J) 

A0(K,J+1)=A0(K,J)+MU*E(K.J)*W(K,J) 

A1(K,J+1)=A1(K,J)+MU*E(K,J)*W(K-1,J) 

A2(K,J+1)=A2(K,J)+MU*E(K,J)*W(K,J-1) 

A3(K,J+1)=A3(K,J)+MU*E(K,J)*W(K-1,J-1) 

50  CONTINUE 

A0(K+1,0)=A0(K,J) 

A1(K+1,0)=A1(K,J) 

A2(K+1,0)=A2(K,J) 

A3(K+1,0)=A3(K,J) 

40  CONTINUE 
C 

C  ******0UTPUT  RESULTS******* 

L=0 

DO  60  K=0,31 
DO  70  J=0,31 
L=L+1 

PRINT  15,  L,E(K,J),A0(K,J),A1(K,J),A2(K,J).A3(K,J) 
WRITE  (80,15)  L,E(K,J),A0(K,J),A1(K,J),A2(K,J),A3(K,J) 
15  FORMAT  ('  ' ,1X,I4,2X,F8.5,2X,F8.5,2X,F8.5,2X,F8.5, 

*  2X,F8. 5) 

70  CONTINUE 

60  CONTINUE 
STOP 
END 
C 

SUBROUTINE  GAUSS(IX,S,AM,V) 

A=0.  0 

DO  50  1=1,12 
CALL  RANDU(IX,IY,Y) 

IX=IY 
50  A=A+Y 

V=(A-6.  0)*S+AM 
RETURN 
END 
C 
C 

SUBROUTINE  RANDU( IX, IY,YFL) 

IY=IX*65539 

IF(IY)5,6,6 

5  IY=IY+2147483647+1 

6  YFL=IY 

YFL=YFL*.  4656613E-9 

RETURN 

END 
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non  o  o  o  n  Cl  n  n  o  rs  n  n  n  n  o  o  o  n  0000 


APPENDIX  B.  LMS  ALGORITHM  IMPLEMENTED  IN  A  NOISE 

CANCELER 

THIS  IS  A  VAX/VMS  FORTRAN  PROGRAM  THAT  IMPLEMENTS  A  IVO 
DIMENSIONAL  2X2  ADAPTIVE  LMS  FILTER  WITHIN  A  NOISE 
CANCELER 

INTEGER  M,N,K,J 

BYTE  B(0: 127) ,BYTEE(0:  127) ,BYTEU(0:  127), 

INTEGER*4  INTE(0: 127,0:  127),INTU(0:  127,0:  127), 

*  IE,IU 

REAL*4  MU,A,FMINE,FMAXE,FMINU,FMAXU 
REAL*4  AA(0:3),E(0:  127,0:  127), U(0:  127,0:  127) 

REAL*4  W(-l: 127,-1: 127), Y(0: 127,0: 127),IM(0: 127,0: 127) 

****^VARIABLE  DEFINITIONS****** 

A=N0ISE  AMPLITUDE 
S=N0ISE  VARIANCE 
AM=N0ISE  MEAN 
AA=FILTER  COEFFICIENTS 
IX=SEED 

MU=SCALING  FACTOR  WITHIN  THE  LMS  ALGORITHM 
IM=INPUT  IMAGE 

W=WHITE  GAUSSIAN  NOISE  GENERATED  BY  SUBROUTINE 
U=SIGNAL  PLUS  N0ISE(IM  +  W) 

Y=FILTER  OUTPUT 
E=ERR0R(U  -  Y) 

FMINE , FMAXE , FMINU , FAMXU=PARAMETERS 

TO  BE  USED  TO  CONVERT  DECIMAL  DATA  TO  BYTE  DATA 

*******inixiaL  VALUES********* 

A=l. 

MU=35.  E-9 
AA(0)=0. 

AA(1)=0. 

AA(2)=0. 

AA(3)=0. 

S=40. 

AM=0. 

IX=65539 
FMINE=1.  E+10 
FMAXE=-1.E+10 
FMINU=1.  E+10 
FMAXU=-1.  E+10 


******0PEN  AN  IMAGE  FILE,  CONVERT  THE  BYTE  DATA  INTO  INTEGERS 
AND  THEN  PLACE  THESE  VALUES  IN  A  MATRIX******* 

OPEN  (UNIT=1,  NAME  ='H0US1G.  DAT' ,  TYPE  ='0LD’,  ACCESS  = 

*  'DIRECT',  RECORDS I ZE=32,  MAXREC=128) 

DO  100  M=0,127 

READ(1'M+1)  (B(J),J=0,127) 

DO  110  N=0.127 
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IF(B(N).LT.  0)THEN 
IM(M,N)=B(N)+256 

PT  ^P 

IM(M,N)=B(N) 

ENDIF 

110  CONTINUE 

100  CONTINUE 
C 

C  *******aDD  white  GAUSSIAN  NOISE  TO  THE  IMAGE  AND  SET  THE 
C  VALUES  OUTSIDE  THE  IMAGE  TO  ZERO******* 

DO  10  M=-l,127 
DO  20  N=- 1,127 

IF  ((M.LT.  0).0R.  (N.LT.O))THEN 
W(M,N)=0. 

IM(M,N)=0. 

ELSE 

CALL  GAUSS(IX,S,AM,V) 

W(M,N)=V*A 

ENDIF 

U(M,N)=IM(M,N)+W(M,N) 

20  CONTINUE 

10  CONTINUE 

C 

C  ******USE  THE  LMS  ADAPTIVE  ALGORITHM  TO  UPDATE 

C  THE  FILTER  COEFICIENTS********** 

DO  40  K=0,127 
DO  50  J=0,127 

Y(K,J)=AA(0)*W(K,J)+AA(1)*W(K-1,J)+ 

*  AA(2)*W(K,J-1)+AA(3)*W(K-1,J-1) 
E(K,J)=U(K,J)-Y(K,J) 

AA( 0 )=AA( 0 )+MU*E( K , J)*W( K , J) 

AA( 1)=AA( 1)+MU*E(K, J)*W(K-1 , J) 
AA(2)=AA(2)+MU*E(K,J)*W(K,J-1) 
AA(3)=AA(3)+MU*E(K,J)*W(K-1,J-1) 

50  CONTINUE 

40  CONTINUE 

C 

C  ********CHANGE  SIGNAL  PLUS  NOISE  AND  ERROR 

C  OUTPUT  INTO  BYTE  DATA  AND  THEN  WRITE 

TO  A  FILE******** 

OPEN  (UNIT=2jNAME=' ERROR.  DAT’ ,TYPE=’NEW' ,ACCESS= 

*  'DIRECT\REC0RDSIZE=32,MAXREC=128) 

OPEN  (UNIT=3,NAME='SIGN0ISE.  DAT’ ,TYPE=’NEW' ,ACCESS= 

*  ’DIRECT' ,REC0RDSIZE=32,MAXREC=128) 

DO  500  1=0,127 

DO  550  J=0,127 

IF(E(I,J).LT.FMINE)THEN 

FMINE=E(I,J) 

ENDIF 

IF(E(  I ,  J).  GT.  FMAXE)THEN 
FMAXE=E(I,J) 

ENDIF 

550  CONTINUE 

500  CONTINUE 

IF  (FMINE. LT. 0)  THEN 
FMAXE=FMAXE-FMINE 
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ENDIF 

DO  600  1=0,127 
DO  650  J=0,127 

IFCFMINE.  LT.  0)THEN 
E(I,J)=E(I,J)-FMINE 
E( I , J)=E( I , J)*255. /FMAXE 
ELSE 

E(I,J)=(E(I,J)=FMINE)*255. /FMAXE 
ENDIF 

1E=NINT(E(I,J)) 

IFCIE.GT. 127)  THEN 
IE=IE-256 
ENDIF 

BYTEE(J)-IE 
650  CONTINUE 

WRITE  (2'I+1)  (BYTEE(N),N=0,127) 

600  CONTINUE 

DO  700  1=0,127 
DO  750  J=0,127 

IFCUC I , J) . LT.  FMINU)THEN 
FMINU=U(I,J) 

ENDIF 

IF(U(  I ,  J).  GT.  FMAXU)THEN 
FMAXU=U(I,J) 

ENDIF 

750  CONTINUE 

700  CONTINUE 

IF  (FMINU. LT.  0)  THEN 
FMAXU=FMAXU-FMINU 
ENDIF 

DO  400  1=0,127 
DO  450  J=0,127 

IF(FMINU.  LT.  0)THEN 
U(I,J)=U(I,J)-FMINU 
U( I , J)=U( I , J)*255.  /FMAXU 
ELSE 

U( I , J)=(U( I , J)=FMINU)*255.  /FMAXU 
ENDIF 

IU=NINT(U(I,J)) 

IF(IU. GT. 127)  THEN 
IU=IU-256 
ENDIF 

BYTEU(J)-IU 
450  CONTINUE 

WRITE  (3'I+1)  (BYTEU(N),N=0,127) 

400  CONTINUE 

CLOSE  (UNIT=2) 

CLOSE  (UNIT=3) 

STOP 

END 

C 

SUBROUTINE  GAUSS(IX,S,AM,V) 

AMP=0. 0 
DO  50  1=1,12 

RANDOM=RAN( IX) 
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AMP=AMP+RANDOM 

CONTINUE 

V=( AMP-6.  0)*S+AM 

RETURN 

END 
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APPENDIX  C.  LMS  ALGORITHM  IMPLEMENTED  IN  AN  ADAPTIVE 

LINE  ENHANCER 

THIS  IS  A  VAX/VMS  FORTRAN  PROGRAM  THAT  IMPLEMENTS  A  TWO 
DIMENSIONAL  2X2  ADAPTIVE  LMS  FILTER  WITHIN  AN  ADAPTIVE 
LINE  ENHANCER 

INTEGER  M,N,K,J 

BYTE  B(0: 127) ,BYTEY(0: 127) ,BYTEU(0: 127) 

INTEGER*4  INTY(0:  127,0:  127),INTU(0:  127,0:  127),IY,IU 
REAL*4  MU,A,FMINY,FMAXY,FMINU,FMAXU 
REAL*4  AA(0:  S^ECO:  127,0:  127), U(0:  127,0:  127) 

REAL*4  W(-l:  127,-1:  127), Y(0:  127,0:  127),IM(0: 127,0:  127) 

***Kr**VARIABLE  DEFINITIONS****** 

A=N0ISE  AMPLITUDE 
S=N0ISE  VARIANCE 
AM=N0ISE  MEAN 
AA=FILTER  COEFFICIENTS 
IX=SEED 

MU=SCALING  FACTOR  WITHIN  THE  LMS  ALGORITHM 
IM=INPUT  IMAGE 

W=WHITE  GAUSSIAN  NOISE  GENERATED  BY  SUBROUTINE 
U=SIGNAL  PLUS  N0ISE(IM  +  W) 

Y=FILTER  OUTPUT 

WW=DELAYED  VERSION  OF  SIGNAL  PLUS  NOISE(U) 

E=ERR0R(U  -  Y) 

FMINY , FMAXY , FMINU , FAMXU=PARAMETERS 

TO  BE  USED  TO  CONVERT  DECIMAL  DATA  TO  BYTE  DATA 


*****Vf* initial  values********* 

A=l. 

MU=35.E-9 

AA(0)=0. 

AA(1)=0. 

AA(2)=0. 

AA(3)=0. 

S=40. 

AM=0. 

IX=65539 
FMINU=1.  E+10 
FMAXU=-1.E+10 
FMINY=1.  E+10 
FMAXY=-1.  E+10 


******0PEN  AN  IMAGE  FILE,  CONVERT  THE  BYTE  DATA  INTO  INTEGERS 
AND  THEN  PLACE  THESE  VALUES  IN  A  MATRIX******* 

OPEN  (UNIT=1,  NAME  ='H0US1G. DAT' ,  TYPE  ='0LD’,  ACCESS  = 

*  'DIRECT' ,  REC0RDSIZE=32,  MAXREC=128) 

DO  100  M=0,127 

READ(1'M+1)  (B(J),J=0,127) 

DO  110  N=0,127 
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IF(B(N).LT.  0)THEN 
IM(M,N)=B(N)+256 

PT 

IM(M,N)=B(N) 

ENDIF 

110  CONTINUE 

100  CONTINUE 
C 

C  *******aDD  white  GAUSSIAN  NOISE  TO  THE  IMAGE  AND  SET  THE 
C  VALUES  OUTSIDE  THE  IMAGE  TO  ZERO******* 

DO  10  M=- 1,127 
DO  20  N=- 1,127 

IF  ((M.LT.  0).OR.  (N.LT.0))THEN 
W(M,N)=0. 

IM(M,N)=0. 

ELSE 

CALL  GAUSS(IX,S,AM,V) 

W(M,N)=V*A 

ENDIF 

U(M,N)=IM(M,N)+W(M,N) 

WW(M.N)=U(M'1,N) 

20  CONTINUE 

10  CONTINUE 
C 

C  **^Hir**usE  THE  LMS  ADAPTIVE  ALGORITHM  TO  UPDATE 

DO  40  K=0,127 
DO  50  J=0,127 

Y(K,J)=AA(0)*WW(K,J)+AA(1)*WW(K-1,J)+ 

*  AA(2)*WW(K,J-1)+AA(3)*WW(K-1,J-1) 
E(K,J)=U(K,J)-y(K,J) 

AA(0)=AA(0)+MU*E(K,J)*WW(K,J) 

AA( 1)=AA( 1)+MU*E(K, J)*WW(K-1, J) 
AA(2)=AA(2)+MU*E(K,J)*WW(K,J-1) 

AA( 3 )=AA( 3 )+MU*E( K , J)*WW(K- 1 , J“ 1 ) 

50  CONTINUE 

40  CONTINUE 

C 

C  ********CHANGE  SIGNAL  PLUS  NOISE  AND  FILTER  OUTPUT 

C  INTO  BYTE  DATA  AND  THEN  WRITE  TO  A  FILE****** 

OPEN  (UNIT=3,NAME=' SIGNOISE.DAT’ ,TYPE='NEW' ,ACCESS= 

*  'DIRECT\REC0RDSIZE=32,MAXREC=128) 

OPEN  (UNIT=4,NAME=' FILTERED. DAT' ,TYPE='NEW' ,ACCESS= 

*  'DIRECT' ,RECORDSIZE=32,MAXREC=128) 

DO  700  1=0,127 

DO  750  J=0,127 

IF(U(I,J).LT. FMINU)THEN 
FMINU=U(I,J) 

ENDIF 

IF( U( I , J) .  GT.  FMAXU)THEN 
FMAXU=U(I,J) 

ENDIF 

750  CONTINUE 

700  CONTINUE 

IF  (FMINU. LT.  0)  THEN 
FMAXU=FMAXU-FMINU 
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ENDIF 

DO  400  1=0,127 
DO  450  J=0,127 

IFCFMINU.  LT.  0)THEN 
U(I,J)=U(I,J)-FMINU 
U( I , J)=U( I , J)*255.  /FMAXU 
ELSE 

U( I , J)=(U( I , J)=FMINU)*255.  /FMAXU 
ENDIF 

IU=NINT(U(I,J)) 

IFCIU.GT. 127)  THEN 
IU=IU-256 
ENDIF 

BYTEU(J)-IU 

CONTINUE 

WRITE  (3'I+1)  (BYTEU(N),N=0,127) 
CONTINUE 
DO  800  1=0,127 
DO  850  J=0,127 

IF( Y( I , J) . LT. FMINY)THEN 
FMINY=Y(I,J) 

ENDIF 

IF(  Y(  I ,  J) .  GT.  FMAXY)THEN 
FMAXY=Y(I,J) 

ENDIF 

CONTINUE 

CONTINUE 

IF  (FMINY.  LT.  0)  THEN 
FMAXY=FMAXY-FMINY 
ENDIF 

DO  900  1=0,127 
DO  950  J=0,127 

IF( FMINY.  LT.  0)THEN 
Y(I,J)=Y(I,J)-FMINY 
Y(I,J)=Y(I,J)*255.  /FMAXY 
ELSE 

Y(I,J)=(Y(I,J)=FMINY)*255.  /FMAXY 
ENDIF 

IY=NINT(Y(I,J)) 

IF(IY. GT. 127)  THEN 
IY=IY-256 
ENDIF 

BYTEY(J)-IY 

CONTINUE 

WRITE  (4'I+1)  (BYTEY(N),N=0,127) 
CONTINUE 
CLOSE  (UNIT=3) 

CLOSE  (UNIT=4) 

STOP 

END 

SUBROUTINE  GAUSS(IX,S,AM,V) 

AMP=0.  0 
DO  50  1=1,12 

RANDOM=RAN( IX) 
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AMP=AMP+RANDOM 

CONTINUE 

V=( AMP-6.  0)*S+AM 

RETURN 

END 
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APPENDIX  D.  RLS  SYSTEM  IDENTIFICATION 

SYSTEM  IDENTIFICATION  OF  AN  ADAPTIVE  RLS  2X2  FILTER  VERSUS 
2X2  RLS  FILTER  VERSUS  A  A  2X2  FILTER  WITH 
KNOWN  COEFFICIENTS 

DIMENSION  W(-l:  32,-1:  32) ,P(0:  3,0:  3),AA(0:  3) ,X(0:  3) 

******VARIABLE  DEFINITIONS****** 

A=N0ISE  AMPLITUDE 
S=N0ISE  VARIANCE 
AM^^NOISE  MEAN 
IX=SEED 

AA=FILTER  COEFFICIENTS 

W=WHITE  GAUSSIAN  NOISE  GENERATED  BY  SUBROUTINE 
D=DESIRED  OUTPUT 
Y=FILTER  OUTPUT 
E=ERR0R(D  -  Y) 

MXM  =  ARRAY  SIZE 
NXN  =  FILTER  SIZE 
Q=WEIGHTING  FACTOR 
P=INVERSE  OF  COVARIANCE  MATRIX 

****^Hf*  INPUT  INITIAL  VALUES  icieirieicirA"k^'h  ick 

M=32 

N=2 

A=l. 

AA(0)=0. 

AA(1)=0. 

AA(2)*0. 

AA(3)=0. 

0=1. 

S=1 

AM=0 

IX=65539 

******  BUILD  INITIAL  ’P'  MATRIX  ********* 

DO  1  K=0,CN**2-1) 

DO  5  J=0,(N**2-1) 

IFCK.EQ.  J)THEN 
P(K,J)=100. 

ELSE 

P(K,J)=0.  0 
ENDIF 
CONTINUE 
CONTINUE 

OPEN  (UNIT=63,FILE='RLSID2D2  DATA* ,STATUS=’ OLD' ) 

*******  CALCULATE  INPUT  MATRIX  (WHITE  GAUSSIAN  NOISE)  ****** 
DO  10  K=-1,M-1 
DO  20  J=-1,M-1 

IF((K.  LT.  0).  OR.  (J. LT. 0))THEN 
W(K,J)=0. 
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ELSE 

CALL  GAUSS(IX,S,AM,V) 

W(K,J)=V*A 

ENDIF 

20  CONTINUE 

10  CONTINUE 

C 

C  *********  CALCULATE  ERROR  BETWEEN  ADAPTIVE  FILTER  OUTPUT 

C  AND  KNOWN  FILTER  OUTPUT.  USE  THIS  ERROR  TO  UPDATE 

C  FILTER  COEFFICIENTS  AND  THE  ’P’  MATRIX  ************* 

DO  40  K=0,M-1 
DO  50  J=0,M-1 
L=K*M+J 

D=.  4*W(K,J)+.  6*W(K-1,J)-.  3*W(K,J-1)+.  2*W(K-1,J-1) 
Y=AA(0)*W(K,J)+AA(1)*W(K-1,J)+AA(2)*W(K,J-1)+ 

*  AA(3)*W(K-1,J-1) 

E=D-Y 

X(0)=W(K,J) 

X(1)=W(K-1.J) 

X(2)=W(K,J-1) 

X(3)=W(K-1,J-1) 

CALL  RLS(P,X,AA,N,E,Q) 

PRINT  15,  L,E,AA(0),AA(1),AA(2),AA(3) 

WRITE  (63,15)  L,E,AA(0),AA(1),AA(2),AA(3) 

50  CONTINUE 

40  CONTINUE 

15  FORMAT  ('  ',2X,I3,2X,F9.5,2X,F9.5,2X,F9.5,2X,F9.5,2X,F9.5) 

STOP 
END 
C 
C 

C  *********  TO  COMPUTE  THE  GAIN  MATRIX  ************ 

SUBROUTINE  RLS(P,X,AA,N,E,Q) 

DIMENSION  P(0:  3,0:  3),X(0:  3),AA(0:  3),TEMY(0: 3,0:  3), 

*  GAMACO:  3),TEMC(0:  3) 

DO  30  L=0,3 

TEM=0. 

DO  20  K=0,3 

20  TEM=TEM+X(K)*P(K,L) 

30  GAMA(L)=TEM 

GAM=Q 

DO  40  K=0,3 

40  GAM=GAM+GAMA(K)*X(K) 

DO  60  L=0,3 
TEM=0. 

DO  50  K=0,3 

50  TEM=TEM+P(L,K)*X(K) 

60  TEMC(L)=TEM 

DO  70  L=0,3 
DO  70  K=0,3 

70  TEMY(L,K)=TEMC(L)*GAMA(K) 

DO  80  K=0,3 
DO  80  L=0,3 

80  P(K,L)=(P(K,L)-(TEMY(K,L)/GAM))/Q 

DO  100  K=0,3 
TEM=0. 
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DO  90  L=Q,3 

90  TEM=TEM+P(K.L)*X(L) 

100  TEMC(K)=TEM 
DO  110  K=0,3 

110  AA(K)=AA(K)+TEMC(K)*E 
Q=.  99*Q+.  01 
RETURN 
END 
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APPENDIX  E.  RLS  ALGORITHM  IMPLEMENTED  IN  A  NOISE 

CANCELER 


THIS  IS  A  VAX/VMS  FORTRAN  PROGRAM  THAT  IMPLEMENTS  A  TVO 
DIMENSIONAL  2X2  ADAPTIVE  RLS  FILTER  WITHIN  A  NOISE 
CANCELER 

INTEGER  M,N,K,J 

BYTE  B(0: 127) ,BYTEE(0:  127) ,BYTEU(0:  127) 

INTEGER*4  INTE(0:  127,0:  127),INTU(0:  127,0: 127),IE,IU 
REAL*4  A,FMINE,FMAXE,FMINU,FMAXU 
REAL*4  AA(0:3),E(0: 127,0:  127), U(0:  127,0:  127) 

REAL*4  P(0:3,0:3),X(0:3),IM(0:127,0:  127)  ( 

REAL*4  W("l: 127,-1:  127), Y(0:  127,0:  127) 

******VARIABLE  DEFINITIONS****** 

A=N0ISE  AMPLITUDE 
S=N0ISE  VARIANCE 
AM=N0ISE  MEAN 
AA=FILTER  COEFFICIENTS 
IX=SEED 

IM=INPUT  IMAGE 

W=WHITE  GAUSSIAN  NOISE  GENERATED  BY  SUBROUTINE 
U=SIGNAL  PLUS  N0ISE(IM  +  W) 

Y=FILTER  OUTPUT 
E=ERR0R(U  -  Y) 

Q=WEIGHTING  FACTOR 
P=INVERSE  OF  COVARIANCE  MATRIX 
FMINE , FMAXE , FMINU , FAMXU=PARAMETERS 

TO  BE  USED  TO  CONVERT  DECIMAL  DATA  TO  BYTE  DATA 
MXM=ARRAY  SIZE 
NXN=FILTER  SIZE 

******nriNITIAL  VALUES********* 

M=128 
N=2 
Q=l. 

A=l. 

AA(0)=0. 

AA(1)=0. 

AA(2)=0. 

AA(3)=0. 

S=40. 

AM=0. 

IX=65539 
FMINE=1.  E+10 
FMAXE=-1.E+10 
FMINU=1.  E+10 
FMAXU=-1.E+10 
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c 

c  ******0PEN  AN  IMAGE  FILE,  CONVERT  THE  BYTE  DATA  INTO  INTEGERS 
C  AND  THEN  PLACE  THESE  VALUES  IN  A  MATRIX******* 

OPEN  (UNIT=1,  NAME  =*H0US1G.  DAT* ,  TYPE  ='OLD',  ACCESS  = 

*  'DIRECT' ,  REC0RDSIZE=32,  MAXREC=128) 

DO  100  K=0,127 

READ(1'K+1)  (B(L),L=0,127) 

DO  110  J=0,127 

IF(B(J).LT.  0)THEN 

IM(K,J)=B(J)+256 

ELSE 

IM(K,J)=B(J) 

ENDIF 

110  CONTINUE 

100  CONTINUE 

CLOSE  (UNITED 
C 

C  ******  build  INITIAL  'P'  MATRIX  ********* 

DO  1  K=0,(N**2-1) 

DO  5  J=0,(N**2-1) 

IF(K.  EQ.  J)THEN 
P(K,J')=100. 

ELSE 

P(K,J)=0.  0 
ENDIF 

5  CONTINUE 

1  CONTINUE 

C 

C  *******adD  WHITE  GAUSSIAN  NOISE  TO  THE  IMAGE  AND  SET  THE 
C  VALUES  OUTSIDE  THE  IMAGE  TO  ZERO******* 

DO  10  K=-1.M-1 
DO  20  J=-1,M-1 

IF  ((K.LT.  0).0R.  (J.LT.  0))THEN 
W(K,J)=0. 

IM(K,J)=0. 

ELSE 

CALL  GAUSS(IX,S,AM,V) 

W(K,J)=V*A 

ENDIF 

U(K,J)=IM(K,J)+W(K,J) 

20  CONTINUE 

10  CONTINUE 

C 

C  **********  CALCULATE  ERROR  BETWEEN  ADAPTIVE  FILTER  OUTPUT 

C  AND  KNOWN  FILTER  OUTPUT.  USE  THIS  ERROR  TO  UPDATE 

C  FILTER  COEFFICIENTS  AND  THE  'P'  MATRIX  ************* 

DO  40  K=0,M-1 
DO  50  J=0,M-1 
L=K*M+J 

Y=AA( 0 )*W( K , J)+AA( 1 )*W( K- 1 , J)+AA( 2 )*W( K , J- 1 )+ 

*  AA(3)*W(K-1,J-1) 

E(K,J)=U(K,J)-Y(K,J) 

EE=E(K,J) 

X(0)=W(K,J) 

X(1)=W(K-1,J) 

X(2)=W(K,J-1) 
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X(3)=W(K-1,J-1) 

CALL  RLS(P,X,AA,N,EE.Q) 

50  CONTINUE 

40  CONTINUE 

********CHANGE  SIGNAL  PLUS  NOISE  AND  ERROR  OUTPUT 
INTO  BYTE  DATA  AND  WRITE  TO  A  FILE***** 

OPEN  (  UNn^2 ,  NAME= '  ERROR.  DAT ' ,  TYPE=  ’  NEW  ’  ,  ACCESS= 

*  'direct' ,REC0RDSIZE=32,MAXREC=128) 

OPEN  ( UNITES, NAME=' SIGNOISE.DAT’ ,TYPE='NEW’ ,ACCESS= 

*  'direct' ,RECORDSIZE=32,MAXREC=128) 

DO  500  1=0,127 

DO  550  J=0,127 

IF(E( I , J). LT. FMINE)THEN 
FMINE=E(I,J) 

ENDIF 

IF(E(  I ,  J).  GT.  FMAXE)'niEN 
FMAXE=E(I,J) 

ENDIF 
550  CONTINUE 

500  CONTINUE 

IF  (FMINE.lt.  0)  THEN 
FMAXE=FMAXE -FMINE 
ENDIF 

DO  600  1=0,127 
DO  650  J=0,127 

IFCFMINE.  LT.  0)THEN 
E(I,J)=E(I,J)-FMINE 
E(I,J)=E(I,J)*255./FMAXE 
ELSE 

E( I , J)=(E( I , J)=FMINE)*255.  /FMAXE 
ENDIF 

IE=NINT(E(I,J)) 

IFCIE.GT. 127)  THEN 
IE=IE-256 
ENDIF 

ByTEE(J)-IE 
650  CONTINUE 

WRITE  (2'I+1)  (BYTEE(N),N=0,127) 

600  CONTINUE 

DO  700  1=0,127 
DO  750  J=0,127 

IF(U( I, J). LT. FMINU)THEN 
FMINU=U(I,J) 

ENDIF 

IF(U(  I ,  J).  GT.  FMAXU)THEN 
FMAXU=U(I,J) 

ENDIF 

750  CONTINUE 

700  CONTINUE 

IF  (FMINU.  LT.  0)  THEN 
FMAXU=FMAXU-FMINU 
ENDIF 

DO  400  1=0,127 
DO  450  J=0,127 

IF( FMINU. LT.  0)THEN 
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U(I,J)=U(I,J)-FMINU 
U( I , J)=U( I , J)*255,  /FMAXU 
ELSE 

U(I,J)=(U(I,J)=FMINU)*255.  /FMAXU 
ENDIF 

IU=NINT(U(I,J)) 

IF(IU,GT.  127)  THEN 
IU=IU-256 
ENDIF 

BYTEU(J)-IU 
450  CONTINUE 

WRITE  (3’I+1)  (BYTEU(N),N=0,127) 

400  CONTINUE 

CLOSE  (UNIT=2) 

CLOSE  (UNIT=3) 

STOP 

END 

C 

c  irHrMrk***  TO  COMPUTE  THE  GAIN  MATRIX  ********** 

SUBROUTINE  RLS(P,X,AA,N,E ,Q) 

DIMENSION  P(0: 3,0:  3),X(0:  3),AA(0:  3),TEMY(0:  3,0:  3), 
*  GAMA(0:3),TEMC(0:3) 

DO  30  L=0,3 
TEM=0. 

DO  20  K=0,3 

20  TEM=TEM+X(K)*P(K,L) 

30  GAMA(L)=^M 

GAM=Q 

DO  40  K=0,3 

40  GAM=GAM+GAMA(K)*X(K) 

DO  60  L=0,3 
TEM=0. 

DO  50  K=0,3 

50  TEM=TEM+P(L,K)*X(K) 

60  TEMC(L)=TEM 

DO  70  L=0,3 
DO  70  K=0,3 

70  TEMY(L,K)=TEMC(L)*GAMA(K) 

DO  80  K=0,3 
DO  80  L=0,3 

80  P(K,L)=(P(K,L''-(TEMy(K,L)/GAM))/Q 

DO  100  K=0,3 
TEM=0. 

DO  90  L=0,3 

90  TEM=TEM+P(K,L)*X(L) 

100  TEMC(K)=TEM 

DO  110  K=0,3 

110  AA(K)=AA(K)+TEMC(K)*E 

Q=. 99*Q+.  01 
RETURN 
END 
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APPENDIX  F.  RLS  ALGORITHM  IMPLEMENTED  IN  AN  ADAPTIVE 

LINE  ENHANCER 


THIS  IS  A  VAX/VMS  FORTRAN  PROGRAM  THAT  IMPLEMENTS  A  TWO 
DIMENSIONAL  2X2  ADAPTIVE  RLS  FILTER  WITHIN  AN  ADAPTIVE 
LINE  ENHANCER 

INTEGER  M.N.K.J 

BYTE  B(0: 127) ,BYTEU(0:  127) ,BYTEY(0:  127) 

INTEGER*4  INTY(0: 127,0:  127),INTU(0:  127,0:  127),IU,IY 

REAL*4  A,FMINU,FMAXU,FMINY,FMAXy 

REAL*4  AA(0:3),E(0:  127,0:  127), U(0: 127,0: 127) 

REAL*4  P(0: 3,0: 3) ,X(0: 3) ,IM(0: 127,0: 127) 

REAL*4  W(-l: 127,-1:  127), Y(0:  127,0:  127) 

******VARIABLE  DEFINITIONS****i«r* 

A=NOISE  AMPLITUDE 
S=N0ISE  VARIANCE 
AM=N0ISE  MEAN 
AA=FILTER  COEFFICIENTS 
IX=SEED 

IM=INPUT  IMAGE 

W=WHITE  GAUSSIAN  NOISE  GENERATED  BY  SUBROUTINE 
U=SIGNAL  PLUS  NOISE(IM  +  W) 

Y=FILTER  OUTPUT 
E=ERROR(U  -  Y) 

WW=DELAYED  VERSION  OF  SIGNAL  PLUS  NOISE(U) 

Q=WEIGHTING  FACTOR 

P=INVARSE  OF  COVARIANCE  MATRIX 

FMINU , FAMXU , FMINY , FMAXY=PARAMETERS 

TO  BE  USED  TO  CONVERT  DECIMAL  DATA  TO  BYTE  DATA 
MXM=ARRAY  SIZE 
NXN=FILTER  SIZE 

******* INITIAL  VALUES********* 

M=128 

N=2 

Q=l. 

A=l. 

AA(0)=0. 

AA(1)=0. 

AA(2)=0. 

AA(3)=0. 

S=40. 

AM=0. 

IX=65539 
FMINU=1.  E+10 
FMAXU=-1.  E+10 
FMINY=1.  E+10 
FMAXY=-1.  E+10 
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c  **^Wr**0PEN  AN  IMAGE  FILE,  CONVERT  THE  BYTE  DATA  INTO  INTEGERS 
C  AND  THEN  PLACE  THESE  VALUES  IN  A  MATRIX******* 

OPEN  (UNIT=1,  NAME  ='H0US1G.  DAT’ ,  TYPE  =’OLD',  ACCESS  = 

*  'DIRECT',  RECORDS I ZE=32,  MAXREC=128) 

DO  100  K=0,127 

READ(1'K+1)  (B(L),L=0,127) 

DO  110  J=0,127 

IF(B(J).LT.  0)THEN 

IM(K,J)=B(J)+256 

ELSE 

IM(K,J)=B(J) 

ENDIF 

110  CONTINUE 

100  CONTINUE 

CLOSE  (UNIT=1) 

C 

c  ******  build  initial  'p'  matrix  ********* 

DO  1  K=0,(N**2-1) 

DO  5  J=0,(N**2-1) 

IFCK.EQ.  J)THEN 
P(K,J)=100. 

ELSE 

P(K,J)=0.  0 
ENDIF 

5  CONTINUE 

1  CONTINUE 

C 

C  *******adD  white  GAUSSIAN  NOISE  TO  THE  IMAGE  AND  SET  THE 
C  VALUES  OUTSIDE  THE  IMAGE  TO  ZERO******* 

DO  10  K=-2,M-1 
DO  20  J=-2,M-1 

IF  ((K.  LT.  0).0R.  (J.  LT.0))THEN 
W(K,J)=0. 

IM(K,J)=0. 

ELSE 

CALL  GAUSS(IX,S,AM,V) 

W(K,J)=V*A 

ENDIF 

U(K,J)=IM(K,J)+W(K,J) 

20  CONTINUE 

10  CONTINUE 

C 

c  *****COMPUTE  THE  DELAYED  VERSION  OF  THE  SIGNAL 

C  PLUS  NOISE  MATRIX(U)****** 

DO  45  K=:- 1,127 
DO  46  J=-l,127 

IF((M.  LT. -l).OR.  (N.LT. -1)  THEN 
WW(K,J)=0. 

ELSE 

WW(K,J)=U(K-1,J) 

ENDIF 

46  CONTINUE 

45  CONTINUE 

C 

C  *********  CALCULATE  ERROR  BETWEEN  ADAPTIVE  FILTER  OUTPUT 
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C  AND  KNOWN  FILTER  OUTPUT.  USE  THIS  ERROR  TO  UPDATE 

C  FILTER  COEFFICIENTS  AND  THE  'P'  MATRIX  **********- 

DO  40  K=0,M-1 
DO  50  J=0,M-1 

Y=AA( 0 )*WW( K . J) +AA( 1 )*WW( K- 1 , J )+AA( 2 ) *WW( K , J- 1 ) + 

*  AA(3)*WW(K-1,J-1) 

E(K,J)=U(K,J)-y(K,J) 

EE=E(K,J) 

X(0)=WW(K,J) 

X(1)^(K-1,J) 

X(2)=WW(K,J-1) 

X(3)=WW(K-1,J-1) 

CALL  RLS(P,X,AA,N,EE,Q) 

50  CONTINUE 

40  CONTINUE 

C 

C  ********CHANGE  SIGNAL  PLUS  NOISE  AND  FILTER  OUTPUT 

C  INTO  BYTE  DATA  AND  WRITE  TO  A  FILE***** 

OPEN  (UNIT=3,NAME=' SIGNOISE.DAT' ,TYPE=’NEW' , ACCESS* 

*  'DIRECT\RECORDSIZE=32,MAXREC=128) 

OPEN  (UNIT=4, NAME*' FILTERED. DAT' .TYPE*' NEW' .ACCESS* 

*  ' direct'. REC0RDSIZE=32.MAXREC*128) 

DO  700  1=0,127 

DO  750  J=0,127 

IF(U(I , J).  LT.  FMINU)THEN 
FMINU=U(I,J) 

ENDIF 

IF( U( I , J) .  GT.  FMAXU)THEN 
FMAXU*U(I,J) 

ENDIF 

750  CONTINUE 

700  CONTINUE 

IF  (FMINU.  LT.  0)  THEN 
FMAXU=FMAXU-FMINU 
ENDIF 

DO  400  1=0,127 
DO  450  J=0,127 

IF(FMINU.  LT.  0)THEN 
U(I,J)=U(I,J)-FMINU 
U(I,J)=U(I,J)*255.  /FMAXU 
ELSE 

U(i,J)=(U(I,J)=FMINU)*255. /FMAXU 
ENDIF 

IU=NINT(U(I,J)) 

IF(IU.  GT.  127)  THEN 
IU=IU-256 
ENDIF 

BYTEU(J)-IU 
450  CONTINUE 

WRITE  (3'I+1)  (BYTEU(N),N=0,127) 

400  CONTINUE 

DO  800  1=0,127 
DO  850  J=0,127 

IF(Y( I , J). LT. FMINY)THEN 
FMINY=Y(I,J) 

ENDIF 
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IF(y(  I ,  J).  GT.  FMAXY)THEN 
FMAXY=Y(I,J) 

ENDIF 

850  CONTINUE 

800  CONTINUE 

IF  (FMINY.  LT. 0)  THEN 
FMAXy=FMAXY-FMINY 
ENDIF 

DO  900  1=0,127 
DO  950  J=0,127 

IFCFMINY.  LT.  0)THEN 
Y(I,J)=Y(I,J)-FMINY 
Y(I,J)=Y(I,J)*255.  /FMAXY 
ELSE 

Y( I , J)=( Y( I , J)=FMINY)*255.  /FMAXY 
ENDIF 

IY=NINT(Y(I,J)) 

IF(IY. GT. 127)  THEN 
IY=IY-256 
ENDIF 

BYTEY(J)-IY 
950  CONTINUE 

WRITE  (4'I+1)  (BYTEY(N),N=0,127) 

900  CONTINUE 

CLOSE  (UNIT=3) 

CLOSE  (UNIT=4) 

STOP 

END 

C 

c  *********  TO  COMPUTE  THE  GAIN  MATRIX  ********** 

SUBROUTINE  RLS(P,X,AA,N,E,Q) 

DIMENSION  P(0;  3,0:  3) .X(0;  3) ,AA(0: 3) ,TEMY(0:  3,0; 3) , 
*  GAMA(0:  3),TEMC(0:  3) 

DO  30  L=0,3 
TEM=0. 

DO  20  K=0,3 

20  TEM=TEM+X(K)*P(K,L) 

30  GAMA(L)=TEM 

GAM=Q 

DO  40  K=0,3 

40  GAM=GAM+GAMA(K)*X(K) 

DO  60  L=0,3 
TEM=0. 

DO  50  K=0,3 

50  TEM=TEM+P(L,K)*X(K) 

60  TEMC(L)=TEM 

DO  70  L=0,3 
DO  70  K=0,3 

70  TEMY(L,K)=TEMC(L)*GAMA(K) 

DO  80  K=0,3 
DO  80  L=0,3 

80  P(K,L)=(P(K,L)-(TEMY(K,L)/GAM))/Q 

DO  100  K=0,3 
TEM=0. 
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DO  90  L=0,3  APP02140 

90  TEM=^M+P(K,L)*X(L)  APP02150 

100  TEMC(K)=TEM  APP02160 

DO  110  K=0,3  APP02170 

110  AA(K)=AA(K)+TEMC(K)*E  APP02180 

Q=.  99*Q+. 01  APP02190 

RETURN  APP02200 

END  APP02210 
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